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Sampling 
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Abstract  The  report  describes  a  Dynamic  Data-Driven  Application  Systems  (DDDAS) 
project  in  which  multiple  mobile  sensors  are  routed  via  a  data-driven  sampling 
scheme.  The  long-term  goal  of  this  project  is  to  provide  a  control-theoretic  frame¬ 
work  to  enable  intelligent,  mobile  systems  to  optimally  collect  sensor-based  obser¬ 
vations  that  yield  accurate  estimates  of  unknown  processes  such  as  aircraft  forma¬ 
tion  flight  and  environmental  monitoring.  The  basic  research  objective  is  to  apply 
tools  from  aerospace  engineering,  specifically  nonlinear  estimation  and  control,  to 
design  coordinated  sampling  trajectories  that  yield  the  most  informative  measure¬ 
ments  of  estimated  dynamical  and  stochastic  systems.  The  technical  approach  to 
achieve  this  objective  is  to  construct  a  framework  for  dynamic,  data-driven  sampling 
algorithms  that  (1)  maximize  the  observability  of  a  nonlinear  dynamical  system 
subject  to  time- varying  perturbations;  and  (2)  minimize  the  uncertainty  in  the  esti¬ 
mate  of  a  nonstationary  random  process  that  requires  nonuniform  sampling.  The  ap¬ 
proach  incorporates  complementary  representations  of  an  unknown  process:  the  first 
uses  a  deterministic,  model-based  parametrization,  whereas  the  second  uses  a  low¬ 
dimensional  statistical  description;  both  approaches  apply  and  enable  the  DDDAS 
concept  in  which  measurement  data  is  used  to  update  the  model  description  and  the 
updated  model  is  used  to  guide  subsequent  measurements. 


1  Introduction 

Observability-based  optimization  in  path  planning  [74]  and  data  assimilation  [38] 
typically  uses  either  a  low-dimensional  parameterized  model  or  an  empirical  data- 
based  representation  of  the  unknown  process;  however,  problems  arise  when  nei¬ 
ther  a  suitable  model  nor  sufficient  data  are  available.  The  novelty  of  the  ap¬ 
proach  described  in  this  report  lies  in  the  study  of  the  observability  of  a  low- 
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Fig.  1  The  DDDAS  framework  for  dynamic,  data-driven  sampling  uses  (left)  information-based 
metrics  to  optimize  sensor  routes,  (middle)  multi-vehicle  control  to  stabilize  the  desired  tra¬ 
jectories,  and  (right)  nonlinear  filtering  to  assimilate  data;  adaptive  sampling  refers  to  the  re¬ 
optimization  of  sensor  routes  that  occurs  after  data  assimilation. 


dimensional  model  of  an  environmental  vector  field  and  the  use  of  the  output  of  a 
data-assimilation  filter  to  guide  the  observability  analysis  via  metrics  from  Bayesian 
experimental  design.  Observability-  and  information-based  optimization  of  sam¬ 
pling  trajectories  has  the  potential  to  yield  a  reliable  and  predictable  capability  for 
intelligent,  mobile  sensors  under  the  Dynamic  Data-Driven  Application  Systems 
(DDDAS)  concept.  Dynamic,  data-driven  Bayesian  nonlinear  filtering  can  exploit 
noisy,  low-fidelity  and  nonlinear  measurements  collected  in  a  distributed  manner 
either  onboard  a  single  sensor  platform  or  by  multiple  sampling  platforms. 

Although  methods  exist  for  the  optimization  of  sampling  trajectories  using  dis¬ 
tributed  parameter  estimation  [17],  optimal  interpolation  [44],  and  heuristic  ap¬ 
proaches  [67],  an  open  question  is  how  to  rigorously  characterize  the  variability 
of  information  content  in  the  unknown  process  and  how  to  target  observations  in 
the  information-rich  regions.  The  merit  of  the  approach  described  in  this  report  lies 
in  the  design  of  a  statistical  framework  based  on  spatiotemporal  estimation  of  non¬ 
stationary  processes  in  meteorology  [35]  and  geostatistics  [33],  and  the  use  of  this 
framework  to  employ  the  DDDAS  concept  for  sampling  environmental  fields.  The 
framework  extends  the  authors’  previous  work  in  this  area  into  multiple  dimensions 
and  incorporates  the  challenge  of  building  a  nonstationary  statistical  representa¬ 
tion  of  a  random  process  while  simultaneously  optimizing  the  sampling  trajectories. 
Figure  1  summarizes  the  overall  workflow  for  data-driven  routing  of  mobile  sensor 
platforms. 

The  power  of  a  DDDAS-based  approach  to  sampling  optimization  is  demon¬ 
strated  by  its  ubiquity  in  fields  ranging  from  robots  for  environmental  monitoring 
and  spatial  estimation  [19,  24,  12,  3,  75,  20]  to  atmospheric  science  and  hurricane 
research  [59,  8,  7,  48,  42].  One  field  that  has  particularly  advanced  the  DDDAS 
framework  is  oceanography  [2,  23,  15,  49,  66,  73,  29,  45,  6],  in  which  autonomous 
underwater  vehicles  are  used  as  mobile  sensors  for  adaptive  sampling.  Indeed  the 
concept  of  optimal  experimental  design  was  first  applied  to  oceanographic  sampling 
in  the  1980’s  [5].  The  optimization  of  sensor  placement  and  data  collection  has  ap- 
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plications  in  fighting  wildfires  [39],  finding  perturbation  sources  in  power  networks 
[1 1],  and  in  spatial  data  collection  for  geostatistics  [52,  70].  Perhaps  it  is  not  surpris¬ 
ing,  given  to  the  range  of  these  applications,  that  there  are  a  variety  of  approaches 
advocated  in  the  literature,  including  adaptation  based  on  maximum  a  posteriori  es¬ 
timation  [72];  stochastic  deployment  policies  [50,  51];  information-based  methods 
[34,  1,  26];  learning  and  artificial  intelligence  [63,  62,  14];  deterministic  methods 
with  heuristic  metrics  [69,  65];  bio-inspired  source  localization  and  gradient  climb¬ 
ing  [25,  27];  and  nonparametric  Bayesian  models  [10].  The  results  described  in  this 
report  differ  from  prior  work  on  adaptive  sampling  of  dynamical  systems  [13]  and 
random  processes  [4,  12],  in  the  novel  application  of  nonlinear  observability  and 
control  coupled  with  recursive  Bayesian  filtering  to  optimize  sensor  routing  for  en¬ 
vironmental  sampling. 

Observability  in  control  theory  is  a  measure  of  how  well  the  state  variables  of  a 
system  can  be  determined  by  measurements  of  its  outputs.  Observability  of  a  linear 
system  is  well  studied  [32]  and  is  characterized  using  the  Kalman  rank  condition, 
which  is  a  special  case  of  the  observability  rank  condition  of  a  general,  nonlinear 
system  [31].  A  nonlinear  system  is  called  observable  if  two  states  are  indistinguish¬ 
able  only  if  the  states  are  identical  [31].  Observability  in  data  assimilation  refers  to 
the  ability  to  determine  the  parameters  of  an  unknown  process  by  using  a  time  his¬ 
tory  of  observations  [37].  Although  standard  observability  tests  give  a  binary  answer 
(observable  or  not),  the  degree  of  observability  can  be  computed  from  the  singular 
values  of  the  observability  gramian  [38],  which  is  a  Hermitian  matrix  of  inner  prod¬ 
ucts  of  the  system’s  outputs.  Computation  of  the  empirical  observability  gramian 
[40]  requires  only  the  ability  to  simulate  the  system,  and  is  therefore  particularly  at¬ 
tractive  for  numerical  optimization.  Algorithms  for  observability-based  estimation 
and  path  planning  have  broader  applications  in  satellite  control  [47]  and  micro  air 
vehicles  [74]. 

Optimal  statistical  interpolation  of  sensor  observations  to  produce  a  stochastic 
estimate  of  an  unknown  random  process  is  based  on  classical  estimation  theory 
[46],  and  is  formerly  known  in  meteorology  and  oceanography  as  objective  analysis 
[9].  Optimal  interpolation  also  yields  a  measure  of  the  uncertainty  or  error  in  the 
estimate,  which  can  be  used  as  a  measure  of  estimator  performance  or  skill  [16]. 
It  is  common  to  compute  estimator  error  under  the  assumption  of  stationarity  of 
the  spatial  and  temporal  variability  of  the  unknown  process  [44],  although  these 
assumptions  may  not  be  borne  out  in  applications  of  interest.  A  stochastic  process 
whose  variability  changes  when  shifted  in  time  or  space  is  called  nonstationary 
[60,  56],  and  methods  exist  to  parametrize  nonstationary  processes  in  oceanography 
[68,  35]  and  geostatistics  [33,  57,  71].  Indeed,  nonstationary-based  strategies  have 
been  previously  applied  to  mobile  sensor  networks  [72],  though  not  based  on  a 
principled  control  design  as  described  here. 

The  remainder  of  the  report  is  organized  as  follows.  Section  2  provides  an 
overview  of  two  complementary  approaches  to  data-driven  adaptive  sampling.  Sec¬ 
tion  3  summarizes  an  application  example  of  observability-based  sensor  routing  for 
aircraft  proximity  flight  as  first  presented  in  [43].  Section  4  reviews  an  application 
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example  of  coverage  control  of  a  nonstationary  spatiotemporal  field  following  [54]. 
Section  5  summarizes  the  report  and  ongoing  research. 


2  Data-driven  Adaptive  Sampling 
2.1  Measures  of  Observability 

Successful  estimation  of  states  in  a  system  depends  on  the  ability  to  observe  the 
desired  states  from  the  sensor  measurements  collected  over  time  [38].  For  instance, 
in  the  application  example  considered  in  Section  3,  one  seeks  to  observe  the  wake 
parameter  states  of  an  aircraft  given  measurements  of  the  differential  pressure  coef¬ 
ficients  along  the  wing  of  a  second  aircraft.  This  section  provides  a  brief  overview 
of  observability  in  the  context  of  linear  and  nonlinear  systems  and  reviews  the  em¬ 
pirical  observability  gramian  [41,  38]  used  to  assess  the  observability  of  the  wake 
parameters  in  the  aerodynamic  model  described  in  Section  3. 

A  dynamical  system  is  said  to  be  observable  if  its  initial  conditions  can  be  de¬ 
termined  from  a  time  history  of  output  measurements  and  control  inputs  over  some 
time  interval  [31,  64].  A  standard  method  of  measuring  observability  is  to  calculate 
the  singular  values  of  the  observability  gramian  [64] .  In  linear  systems  theory,  the 
singular  values  of  the  observability  gramian  quantify  the  ease  in  determining  the 
initial  states  from  the  outputs  generated  over  time  as  follows:  large  singular  values 
imply  the  mapping  is  easily  invertible,  whereas  small  or  zero  singular  values  imply 
it  is  not  [64,  p.  125-126].  This  report  quantifies  the  unobservability  of  a  system  using 
the  unobservability  index  ,  which  is  the  reciprocal  of  the  smallest  singular  value 
Gmin->  he.,  [38] 

«  =  0) 

The  unobservability  index  (1)  reflects  the  least  observable  mode  in  the  system  and 
provides  a  worst-case  observability  measure  for  the  system. 

The  observability  of  a  nonlinear  system  may  be  difficult  to  determine  analyt¬ 
ically,  because  it  requires  tools  from  differential  geometry  [31].  If  the  dynamical 
model  of  interest  is  solved  numerically,  it  is  justified  to  pursue  numerical  tech¬ 
niques  for  calculating  the  nonlinear,  empirical  observability  gramian  [38,  28].  The 
empirical  observability  gramian  does  not  require  linearization,  which  may  fail  to 
adequately  model  the  input/output  relationship  of  the  nonlinear  system  over  a  wide 
range  of  operating  conditions,  but  merely  the  ability  to  simulate  the  system  [38]. 
Indeed,  the  empirical  observability  gramian  maps  the  input-output  behavior  of  a 
nonlinear  system  more  accurately  than  the  observability  gramian  produced  by  lin¬ 
earization  of  the  nonlinear  system  [28].  The  empirical  observability  gramian  is  a 
square  matrix  whose  dimension  matches  the  size  of  the  state  vector  and  whose  /,  jth 
entry  represent  the  sensitivity  of  the  output  to  infinitesimal  perturbations  about  their 
nominal  value  of  the  corresponding  i  and  j  states  or  unknown  parameters. 
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Fig.  2  The  log  of  the  unobservability  index  of  an  aircraft  wake  over  a  range  of  relative  positions  of 
a  follower  aircraft  equipped  with  distributed  pressure  sensors.  Reprinted  from  [43]  with  permission 
from  the  American  Institute  of  Aeronautics  and  Astronautics. 

The  empirical  observability  gramian  is  formally  defined  as  follows.  Let  £,-e/, 
i  =  1 , . . . ,  M,  be  a  small  displacement  of  the  nominal  parameter  along  the  ith  unit 
vector  e*  G  and  let  Q  G  Mm  be  the  set  of  nominal  parameter  values.  The  (/,  y)th 
component  of  the  M  x  M  empirical  observability  gramian  Wo  is  [38] 


(2) 


i  1 , . . . ,  M,  j  1 , . . . ,  Af , 


where  Q±l  =  Q  ±  £;-e/  produces  the  output  =  h(,0±i,u).  The  observability  of 
a  nonlinear  system  is  measured  by  calculating  the  unobservability  index  (1)  of  the 
empirical  observability  gramian  Wo . 

Figure  2  depicts  the  observability  of  an  aircraft  wake  in  a  plane  orthogonal  to  a 
follower  aircraft,  revealing  blind  spots  in  which  the  follower  aircraft  may  not  be  able 
to  detect  the  leader’s  wake.  Evaluation  of  the  empirical  observablity  gramian  over 
the  space  of  leader  aircraft  positions  provides  an  observability  map  of  the  leader 
aircraft  wake  (see  Section  3  for  model  details).  The  observability  analysis  is  used  in 
Section  3  to  guide  the  design  of  a  recursive  Bayesian  filter  for  estimating  the  leader’s 
wake  parameters  and  an  observability -based  control  strategy  to  steer  the  follower  to 
a  desired  relative  position.  Maintaining  adequate  observability  along  the  trajectory 
guarantees  the  performance  of  the  filtering  scheme. 


2.2  Evaluating  Mapping  Error 

A  spatiotemporal  field  is  statistically  described  by  its  mean  and  the  covariance  func¬ 
tion  between  any  two  points  i  and  j.  A  covariance  function  is  a  positive-definite 
function  that  describes  the  variability  of  the  field  between  the  ith  and  yth  loca¬ 
tion,  as  described  in  [6].  A  field  is  stationary  if  its  covariance  function  depends 
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Fig.  3  Two  vehicles  traveling  through  (a)  stationary  and  (b)  nonstationary  spatiotemporal  fields. 
The  tubes  represent  the  volume  covered  by  sensor  measurements.  Reprinted  from  [54]  with  per¬ 
mission  from  Elsevier. 


only  on  the  relative  position  of  i  and  j  and  and  is  nonstationary  if  it  depends  on  i 
and  j  independently.  There  are  a  number  of  choices  for  the  form  of  a  nonstation¬ 
ary  covariance  function,  e.g.,  Matern,  rational  quadratic,  Ornstein-Uhlenbeck  and 
squared-exponential  forms  as  described  in  [33].  For  the  statistics-based  sampling 
strategy  described  in  Section  4,  we  require  a  covariance  function  that  is  a  product  of 
a  spatial  covariance  function  and  a  temporal  covariance  function. 

To  represent  spatiotemporal  fields  with  non-uniform  coverage  requirements,  we 
adopt  a  nonstationary  squared  exponential  covariance  function  of  the  following 
form,  as  introduced  in  [33]: 


C(ri,rj ) 


|£(n)|1/4|£(o)l1/4 


Z(.n)+Z(rj) 

2 


1/2 


exp 


(3) 


where  Z(r^)  G  M3x3  is  a  positive  definite  symmetric  matrix  that  is  continuous  in 
the  position  r *  G  M2.  The  square  roots  of  the  diagonal  elements  of  Z(t>)  are  the 
spatial  and  temporal  decorrelation  scales  of  the  field.  The  decorrelation  scales  are 
the  spatial  and  temporal  separations  at  which  the  covariance  function  evaluates  to 
l/e  «  0.368.  Note,  for  a  stationary  field,  the  decorrelation  scales  are  constant,  but 
for  a  nonstationary  field  they  may  vary  in  space  and  time.  The  covariance  func¬ 
tion  (3)  is  used  to  derive  the  coordinate  transformation  in  Section  4,  which  clusters 
measurements  in  space-time  regions  with  short  decorrelation  scales,  and  spreads 
measurements  elsewhere. 


FA9550-1 3-1-0162 


Data-driven  Routing  of  Autonomous  Vehicles 


7 


The  statistics-based  sensor-routing  application  of  the  DDDAS  concept  seeks  to 
provide  optimal  coverage  of  an  estimated  spatiotemporal  field.  The  coverage  is 
deemed  optimal  when  the  measurement  density  in  space  and  time  is  proportional 
to  the  variability  of  the  field.  To  determine  when  measurements  are  redundant,  we 
consider  the  footprint  of  a  measurement,  defined  as  the  volume  in  space  and  time 
contained  in  an  ellipsoid  centered  at  the  measurement  location  with  principle  axes 
equal  to  the  decorrelation  scales  of  the  field.  Figure  3  depicts  two  vehicles  taking 
measurements  along  a  circular  trajectory  in  the  space-time  domain.  The  red  and  blue 
tubes  along  the  vehicles’  trajectories  are  the  sensor  swaths  created  by  the  measure¬ 
ment  footprints.  Figure  3(a)  shows  the  swaths  for  a  stationary  field  and  Figure  3(b) 
shows  the  swaths  for  a  nonstationary  field  in  which  the  temporal  scales  contract  at 
t  =  10  hours.  The  goal  is  to  design  the  vehicle  trajectories  so  that  the  swaths  created 
by  the  set  of  all  measurement  footprints  cover  the  entire  field  with  minimal  overlaps 
or  gaps,  even  when  the  decorrelation  scales  of  the  field  vary. 

To  determine  the  mapping  error,  we  employ  optimal  interpolation  as  in  [6] .  Con¬ 
sider  the  field  A(r)  to  be  a  discrete  random  field  in  space  r  =  (x,y)  and  time  t.  Let 
fd  =  [xd,  yd,  td]T  be  the  space-time  location  of  measurement  d  =  1  where 

D  is  the  total  number  of  measurements  taken  by  all  of  the  vehicles,  and  £d  be  the 
measurement  noise,  so  that  the  value  of  measurement  d  is  Zd  =  A{rd)  +  £d.  We  as¬ 
sume  E[£m£j]  =  do Smi,  where  E[]  denotes  the  expected  value,  do  is  the  standard 
deviation  of  the  measurement  noise,  and  <5m/  is  the  Kronecker  delta,  which  implies 
that  the  noise  from  any  two  distinct  measurements  is  uncorrelated  and  the  variance 
of  the  noise  is  do. 

Let  Ce  denote  the  covariance  of  the  error  in  the  estimate  A(r)  after  assimilating 
the  set  of  measurements  f  =  [r\ , . . . ,  ro]T •  We  have 

D  D 

Ce{ri,rj;r)=C(ri,rj)-  £  ^C^ridd)(M~l)dlC(rj,n), 
d=\l=\ 

as  explained  in  [22,  9],  where  M~l  is  the  inverse  of  the  measurement  covariance 
matrix  whose  entries  are  [44] 


Mdi=E[zdZi\  =C(rd,ri)  +  d0Sdi. 


The  mapping  error  is  the  diagonal  of  the  error  covariance  matrix  Ce.  The  aver¬ 
age  (resp.  maximum)  mapping  error  is  computed  by  averaging  (resp.  finding  the 
maximum  of)  all  of  the  elements  of  the  mapping  error.  Since  the  vehicles  sample 
uniformly  in  time,  the  mapping  error  is  minimized  in  a  stationary  field  by  traveling 
at  maximum  speed  to  place  as  many  measurements  as  possible  in  the  domain,  as 
further  illustrated  in  [58,  53].  Traveling  as  fast  as  possible  may  not  be  the  best  solu¬ 
tion,  however,  for  a  nonstationary  field,  since  one  needs  to  slow  down  to  concentrate 
more  measurements  in  a  space/time  region  of  high  variability. 
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3  Application  I:  Wake  Estimation  and  Formation  Control 


3.1  Aerodynamic  Model 


This  section  illustrates  the  DDDAS  concept  using  an  observability-based  sensor 
routing  application  originally  presented  in  [43].  Consider  two  aircraft  in  steady  level 
flight  through  an  inviscid,  incompressible,  irrotational  fluid.  We  employ  an  aerody¬ 
namic  model  of  one  aircraft  flying  in  the  other  aircraft’s  wake,  similar  to  that  of 
Hemati  et  al.  [30]  and  Pachter  et  al.  [55].  Let  the  reference  frame  =  (O,  &i ,  &2,  £>3) 
with  origin  O  be  attached  to  the  center  of  the  leading  edge  of  the  follower’s  wing 
with  basis  vectors.  Assume  that  the  follower  aircraft  maintains  kinematic  control 
of  its  vertical  and  horizontal  velocities  such  that  the  velocity  of  frame  with  re¬ 
spect  to  the  lead  aircraft  in  steady  level  flight  is  V/  =  Vyb2  +  Vzbi  (the  bi  com¬ 
ponent  is  assumed  to  be  zero).  The  leader  has  wingspan  b  and  center  position 
Vl  =  XLb\  +  yi)>2  +  ZLbi  relative  to  O.  Assume  \xl\  is  sufficiently  large  (i.e.,  greater 
than  two  wingspans  [55])  such  that  the  wake  of  the  leader  is  adequately  represented 
using  potential  flow  theory  as  the  sum  of  two  infinite-line  vortices,  each  with  cir¬ 
culation  strength  fl,  extending  horizontally  behind  the  wingtips  of  the  leader  along 
the  b\  direction.  The  Biot-Savart  law  gives  the  following  vertical  component  of  the 
wake  [55,  30]  wl  at  a  point  (x,y,z)  along  the  b2  axis  (x  =  z  =  0)  as  a  function  of  the 
leader’s  position1 : 


wL(y,rL,yL,ZL,b ) 


rL(y-yL-b/2) _ rL{y-yL+b/2) 

27i(zl+(y-yL-b/2)2)  27i{zl+(y-yL+b/2)2) ' 


(4) 


Note  that  (4)  is  symmetric  about  zl  =  0-  This  property  necessitates  invoking  a  sec¬ 
ond  sensor  to  break  the  vertical  symmetry.  The  aerodynamic  signature  on  the  fol¬ 
lower  aircraft  created  by  the  upwash  field  of  the  leader  will  be  used  to  estimate  the 
two-dimensional  position  r  =  of  the  leader  relative  to  the  follower  and  the 

circulation  strength  Tl  of  the  leader’s  wake.  The  position  estimate  r  will  be  used  in 
an  optimal  controller  to  steer  the  follower  to  the  desired  relative  position. 

The  flow  around  the  follower  aircraft  in  response  to  the  upwash  field  of  the  leader 
is  modeled  by  employing  the  lifting-line  solution  following  Katz  and  Plotkin  [36, 
p.331-340].  For  simplicity,  assume  the  follower  is  represented  by  a  flat,  thin,  rect¬ 
angular  wing  with  large  aspect  ratio  ( M  >  4)  and  chord  length  c.  Since  the  fluid 
is  inviscid,  incompressible,  and  irrotational,  its  motion  can  be  represented  by  the 
gradient  of  a  potential  function  that  satisfies  Laplace’s  equation  [36] 


V20  =  O.  (5) 

In  addition,  at  every  point  on  the  wing  the  potential  function  <P  must  satisfy  a 
boundary- value  constraint  that  ensures  that  there  is  no  flow  normal  to  the  wing  sur¬ 
face.  Assuming  the  wing  is  thin,  its  normal  vector  n;  at  any  point  (x^yt^Zi)  on  the 


1  The  notation  g(a,b\a,ft)  represents  a  function  g(-)  that  depends  on  the  state  variables  a,b  and 
the  parameters  a  and  /3 . 
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Fig.  4  Leader  and  follower  aircraft  representations.  The  wake  of  the  leader  produces  an  aerody¬ 
namic  signature  on  the  follower  through  its  upwash  field  (blue).  Reprinted  from  [43]  with  permis¬ 
sion  from  the  American  Institute  of  Aeronautics  and  Astronautics. 


wing  surface  is  approximately  n,-  «  b3,  i  =  1, . . .  ,N,  which  implies 

V0b3  =  O.  (6) 

To  satisfy  these  constraints,  lifting-line  theory  constructs  a  suitable  potential  func¬ 
tion  from  a  collection  of  line  vortices.  N  equally  spaced  horseshoe  vortices  are 
bound  to  the  quarter  chord  of  the  follower  aircraft  wing  such  that  2N  trailing  vor¬ 
tices  extending  infinitely  downstream.  The  zth  bound  horseshoe  vortex  has  circula¬ 
tion  strength  J].  The  number  N  of  horseshoe  vortices  must  be  chosen  large  enough 
for  adequate  model  fidelity,  yet  small  enough  to  remain  computationally  tractable. 
The  freestream  fluid  velocity  Uoo  has  magnitude  £/«,  and  angle  of  attack  a  relative 
to  the  wing.  (Assume  the  freestream  velocity  has  zero  sideslip  and  a  is  small.) 

Since  the  line  vortex  is  a  solution  to  Laplace’s  equation  (5)  [36],  the  flow  due  to 
the  freestream  velocity,  leader  aircraft  upwash,  and  horseshoe  vortices  must  satisfy 
the  normal  flow  constraint  (6).  Therefore,  (6)  evaluated  at  any  given  point  on  the 
wing  must  satisfy 

Whs  +  wl  +  Vz  + Uoo  sin  cl  —  0,  (7) 

where  Whs  is  the  b3  component  of  the  flow  generated  by  the  horseshoe  vortices,  wl 
is  the  vertical  component  of  the  leader  aircraft’s  wake  given  by  (4),  Vz  is  the  b3  com¬ 
ponent  of  the  inertial  velocity  of  the  follower  aircraft  expressed  in  frame  and  the 
fourth  term  on  the  left-hand  side  is  the  normal  component  of  the  freestream  veloc¬ 
ity.  Note  that  this  model  neglects  aerodynamic  influences  due  to  aircraft  pitching, 
rolling,  and  sideslip  maneuvers,  under  the  assumption  that  in  close  proximity  these 
motions  are  negligible  relative  to  the  aerodynamic  effect  of  vertical  motion. 


FA9550-1 3-1-0162 


10 


Derek  A.  Paley 


To  solve  for  the  horseshoe  vortex  strengths  I]  that  satisfy  (7),  we  employ  the  col¬ 
location  method  following  Katz  and  Plotkin  [36,  p.331-334].  Impose  the  constraint 
(7)  at  N  collocation  points  centered  at  each  horseshoe  vortex  along  the  3 /4-chord 
line.  Since  the  flow  at  any  collocation  point  is  linearly  dependent  on  the  circulation 
strength  J]  of  the  ith  horseshoe  vortex,  (7)  applied  at  the  N  collocation  points  forms 
a  set  of  N  linear  algebraic  equations  with  N  unknown  circulation  strengths  [36]  J], 

an  •••  a\M  7] 

aNi  •••  aNN J  |/w 

Here 

N 

Y,a]iri  =  whs{yj ),  (9) 

i=  1 

is  the  normal  component  of  the  flow  at  the  j  collocation  point,  Vz  is  the  (control) 
velocity  of  the  follower  aircraft,  and  WL{yj)  is  obtained  by  evaluating  (4)  at  yj. 
Equation  (8)  is  solved  for  the  circulation  strength  distribution  I"j , . . . ,  by  inverting 
the  aji  coefficient  matrix  [36]. 

The  circulation  strength  distribution  r  (y )  =  lim^^co  I]  is  used  to  calculate  mea¬ 
surable  quantities  that  will  be  used  in  the  nonlinear  wake  estimation  process.  For 
example,  Hemati  et  al.  [30]  use  the  following  measurements  of  the  differential  pres¬ 
sure  coefficient  ACp\ 

ACP(x,y)==^(l- 1)1/2.  (10) 

Assume  the  measurements  of  differential  pressure  are  equally  spaced  on  the  3 74- 
chord  line  v  =  (3c) /4  along  the  span  of  the  wing. 

Let  h (Cl)  =  [ACp\ . . .  ACpp]T  beaPx  1  column  matrix  of  differential  pressure 
measurements  calculated  using  (10)  on  the  3 /4-chord  line,  where  P  is  the  number  of 
measurements.  Taking  h(Cl)  as  the  output  and  using  the  leader  aircraft  dynamics  in 
frame  generates  the  following  state-space  form  of  the  input-output  relationship 
between  the  wake  parameter  states  Cl  =  \yL,ZL,I V\T  and  the  measurements  (j): 


1 

wl(ji) 

=  —  (f/oosina  +  VZ) 

1 

— 

wl(jn) 

~-Vy~ 

n=  -vz 

o 

0  =  h(fl). 


(11) 


The  model  (11)  is  used  to  evaluate  the  observability  of  the  state  Cl  with  the  output 
equation  h(fl),  design  an  observer  to  estimate  Cl  from  noisy  output  measurements, 
and  implement  an  optimal  control  strategy  for  the  follower  aircraft. 
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Fig.  5  A  schematic  diagram  of  an  observability-based  sampling  algorithm.  A  recursive  Bayesian 
filter  provides  parameter  estimates  Q.  from  noisy  measurements  a.  The  estimated  parameters  are 
used  to  calculate  observability  optimizing  control  parameters  % *  that  characterize  the  multi-sensor 
sampling  formation. 


3.2  Formation  Flight  Control 

The  optimal  control  approach  described  here  incorporates  a  weighted  wavefront 
expansion,  known  as  the  fast  marching  method  [61],  to  generate  an  optimal  cost-to- 
go  potential  map  relative  to  the  desired  position  of  the  lead  aircraft.  The  gradient  of 
this  potential  function  provides  the  optimal  control  with  respect  to  an  observability- 
based  cost  function.  This  method  does  not  require  iterative  calculation  of  the  control 
since  the  potential  is  calculated  over  the  entire  relative  position  space,  making  the 
control  computationally  less  expensive  than  a  receding-horizon  control  [43]. 

The  optimal  path  planning  problem  is  formulated  as  follows  [61,  pg.  284-291]. 
Given  a  desired  position  r des,  the  goal  is  to  find  the  path  J £?(/):  [0,°o)  -A  M2  from 
Ydes  to  any  point  ro  that  minimizes  the  observability -based  cost  integral  [61] 

f  log^(/))d/, 

^ rdes 

where  /  is  the  arc-length  parameterization  of  the  path  ££  and  £  (•)  is  the  unobserv¬ 
ability  index  evaluated  along  ££ .  Let  the  minimum  cost  required  to  travel  from  r des 
to  a  point  r  be  [61] 

Jwf( r)=min  f  §(jg?(Z))d/,  (12) 

^  Y  des 

such  that  the  level  set  Jwf{?)  =  C  is  the  set  of  points  that  can  be  reached  with 
minimal  cost  C.  By  construction,  level  sets  are  orthogonal  to  the  minimal  cost  paths 
[61]  implying  that  the  optimal  path  descends  the  gradient  of 

The  fast  marching  method  [61]  is  a  wavefront  propagation  technique  that  is  used 
to  efficiently  compute  Jwf(t)  for  the  domain  around  the  leader  aircraft  [61,  pg.  86- 
99].  Since  the  optimal  path  descends  the  gradient  of  7^(r),  an  observer-based  feed- 
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(a) 


Fig.  6  Simulation  illustrating  aerial  refuel  positioning  with  optimal  control  algorithm  using 
wavefront-propagation  method,  (a)  Trajectory;  (b)  marginal  probability  densities  of  Bayesian  filter. 
Reprinted  from  [43]  with  permission  from  the  American  Institute  of  Aeronautics  and  Astronautics. 


5  10 

Time  (s) 


5  10  15 

Time  (s) 

(C) 


back  control  incorporates  estimates  of  the  relative  position  r  according  to 

u  =  (13) 

where  the  gain  Kwf  >  0.  Assuming  the  desired  final  location  is  fixed,  the  cost  po¬ 
tential  JW(r)  need  be  calculated  only  once  to  produce  all  possible  optimal  paths. 

Numerical  simulations  of  the  wavefront-propagation  control  strategy  were  com¬ 
puted  for  the  follower  aircraft  positioning  itself  for  aerial  refueling.  The  simula¬ 
tion  assumes  noise  E  =  10-5  in  the  differential  pressure  coefficient  measurements 
and  total  simulation  time  T  =  15  seconds.  A  recursive  Bayesian  filter  assimilates 
measurements  at  5  Hz  [43].  The  wavefront-propagation  control  strategy  has  gain 
Kwf  =  5  for  all  simulations.  In  aerial  refueling  the  follower  aircraft  intercepts  a 
filling  nozzle  that  extends  outward  from  the  tail  of  the  leader  [18]  by  maintaining 
a  desired  relative  position  of  r des  —  (0,0)  wingspans.  This  implementation  assumes 
the  aerodynamic  effects  of  the  filling  nozzle  on  the  follower  aircraft  are  minimal  at 
the  sensor  positions.  Figure  6  depicts  the  simulations  results,  which  show  the  trajec¬ 
tory  converging  to  the  desired  position  along  a  path  of  guaranteed  observability  of 
the  leader  aircraft’s  relative  position. 


4  Application  II:  Coverage  Control  of  a  Nonstationary  Field 
4.1  Coordinate  Transformation 

This  section  defines  a  nonlinear  coordinate  transformation  that  renders  a  nonstation¬ 
ary  field  locally  stationary,  following  [54].  The  significance  of  this  transformation 
is  that  it  permits  the  use  of  multi- vehicle  sampling  algorithms  designed  for  use  in 
stationary  fields.  We  derive  two  transformations,  one  for  which  the  spatial  decorre- 
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lation  scales  are  coupled,  i.e.,  ox  =  ox(xk,yk)  and  oy  =  ay  (xk,yk),  and  another  for 
which  the  spatial  scales  are  decoupled,  i.e.,  ax  =  <Jx(xk)  and  ay  =  oy(yk). 

Let  rk  =  [xk,  yk,  tk]T  be  the  space-time  coordinates  in  the  original  domain,  called 
the  r-domain,  and  Rk  =  [Xk,  Yk,  Tk]r  be  a  set  of  transformed  coordinates  in  a  new 
domain,  called  the  7?-domain.  Also,  let  G(V,E)  be  a  lattice  graph  with  nodes  V 
placed  at  each  Rk  and  undirected  edges  E  connecting  adjacent  cells  (the  four  nearest 
neighbors  in  two  dimensions).  The  lattice  graph  ensures  that  adjacent  nodes  are 
locally  stationary  under  the  transformation.  We  would  like  to  find  Rk  such  that  the 
field  with  covariance  (3)  becomes  locally  stationary,  i.e.,  all  adjacent  nodes  in  G 
are  (approximately)  stationary,  as  in  [60].  Intuitively,  this  corresponds  to  ensuring 
that  the  field  is  stationary  in  the  neighborhood  of  each  point  Rk.  Let  £(r;-,r;)  = 
(E(ri)+E(rj))/ 2,  r*/  =  n  —  rj ,  and  Rij  =Ri~Rj.  Thus,  we  seek  to  satisfy 

^(^)|(/4^(0)|1/4J—  rjj^))-^]  =e[-RjjRiJ]  (14) 

y/\Z(n,rj)\ 

between  all  adjacent  points  i  and  j.  The  left-hand  side  of  (14)  is  the  nonstation¬ 
ary  covariance  function  in  the  r-domain  and  the  right  hand  side  is  the  covariance 
function  in  the  7?-domain,  which  we  desire  to  be  stationary.  Without  loss  of  gener¬ 
ality,  the  spatial  and  temporal  decrorrelation  scales  in  the  ^-domain  are  chosen  to 
be  Ox  =  oy  =  Tt  =  1 .  Since  the  spatial  decorrelation  scales  are  coupled,  the  trans¬ 
formations  are  also  coupled;  however,  the  temporal  component  of  (14)  is  decoupled 
by  assumption.  Let  axi  =  ex*  (*,-,?,•),  ayi  =  oy(xi,yi ),  and  t,-  =  T (f,-).  Then,  (14)  yields 
the  following  two  relations: 


(Xi-Xj)2  +  (Yi-Yj)2  = 


+  ln 


(: Ti-Tj y  = 


t?  +  t] 


(15) 

(16) 


Equations  (15)  and  (16)  show  how  the  space-time  separation  between  two  points 
in  the  ^-domain  depends  on  the  nonstationarity  in  the  r-domain.  However,  to  com¬ 
plete  the  spatial  transformation  we  need  another  relationship  between  X(j  =  Xj  —  Xj 
and  Yij  =  Yt-  Yj  to  separate  (15)  into  two  equations.  We  use  the  orientation  between 
points  in  the  7?-domain  and  r-domain.  The  relation 


tan 


Yi-Yj 

Xi-X: 


=  tan 


yi-yj 


(17) 


preserves  the  relative  orientation  between  points  i  and  j ,  which  ensures  that  the 
transformation  will  not  rotate  the  coordinates  in  the  ^-domain,  as  it  does  in  [60]. 
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Substituting  (17)  in  (15)  we  obtain 

2  (xj-xj)2  2  (yi-yj)2 


Xu  =  ± 


Yij  =  ± 


°xi  +  °xj 


■  +  ■ 


ayi  +  <yyJ 


+  ln 


yt-yj 


-  dx,ij, 


2  (xi-Xj)2  2  (yi-yj)2 


(J.2  +  <7? 

i  1  J 


CJ-2  +  CJ? 

'  1  J 


+  ln 


i  +  _ _ xi- 

yt-yj 


—  dY,ij ,  and 


Ti-Tj  =  ± 


2  (A  Q')2 

t2  +  t2 


-In 


2ZiZj 


2oxjGxj 


2oxioXj 


1/2 


~  dT,ij- 


iGyiGyj 


2  Gyi  Gy  j 


1/2 


(18) 


1/2 


(19) 

(20) 


Equations  (18)-(20)  represent  the  transformation  from  the  r-domain  to  the  /?- 
domain.  Note  that  the  choice  of  plus  or  minus  in  (18)-(20)  is  arbitrary  as  it  rep¬ 
resents  flipping  the  new  coordinates  about  the  corresponding  axis.  Thus,  we  use  the 
positive  root. 

Let  Z?  be  the  incidence  matrix  of  the  lattice  graph  G.  We  rewrite  (18)-(20)  as 

BT X  =  dx ,  BTY  =  dY,  and  BrT  =  dT ,  (21) 


where  dx  is  the  vector  of  all  terms  dxjj  and  dy  and  dj  are  defined  similarly.  These 
equations  are,  in  general,  overdetermined.  The  least-squares  solution  to  (21)  is  ob¬ 
tained  using  the  pseudoinverse.  The  use  of  the  lattice  graph  topology  is  justified 
since  points  with  short  separation  (in  space  or  time)  have  a  stronger  impact  on 
the  mapping  error  than  points  with  large  separation.  We  have  also  found  the  lat¬ 
tice  topology  to  be  particularly  sensitive  to  nonstationarities  in  the  r-domain;  higher 
connected  graphs  tend  to  concentrate  measurements  less. 

The  inverse  transformation  from  the  ZGdomain  to  the  r-domain  is  solved  by  cre¬ 
ating  a  lookup  table  as  follows.  A  uniform  grid  is  created  in  the  r-domain  and 
mapped  to  the  ZGdomain,  which  produces  a  nonuniform  grid  of  locations  in  the 
ZGdomain.  Trajectories  generated  in  the  ZGdomain  are  discretized  and  mapped  back 
to  the  r-domain  by  interpolating  this  grid.  Figure  7  shows  an  example  of  the  two- 
dimensional  spatial  transformation  using  decorrelation  scales  crx  and  CJy  and  a  vary¬ 
ing  spatial  decorrelation  scale  where  a  “dip”  in  the  scales  is  centered  at  (3,3).  On 
the  right  is  the  r-domain,  which  has  a  square  boundary.  When  the  r-domain  bound¬ 
ary  is  mapped  to  the  ZGdomain,  it  generates  a  new,  curved  boundary.  A  uniform  grid 
mapped  from  the  ZGdomain  to  the  r-domain  shows  how  space  is  warped  to  concen¬ 
trate  measurements  near  the  nonstationarity.  For  the  conditions  that  guarantee  the 
mapping  is  invertible,  the  reader  is  referred  to  [54]. 

Next,  consider  when  the  spatial  decorrelation  scales  are  decoupled,  i.e.,  crx?  = 
Gx(xi)  and  Gyi  =  Gy(yi).  In  this  case,  (18)-(20)  reduce  to 
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-5  0  5 
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Fig.  7  Example  of  the  spatial  transformation  with  a  nonstationarity  at  (3,3).  Reprinted  from  [54] 
with  permission  from  Elsevier. 


Z-Xj  = 


Yi~Yj  = 


Ti  —  Tj  = 


2(x,-—  xjf 

2  (yi-yj)2 


+  ln 


+  ln 


yi  '  yj 

Srf  +  ln 


,  and 


(22) 

(23) 

(24) 


For  a  lattice  topology  with  decoupled  decorrelation  scales,  we  have  j  =  i  +  1,  i.e., 
Xj  =  Xi  +  AXi9  xj  =  Xi  +  (and  similarly  for  7,  y,  T,  and  0-  Taking  the  limit 
as  the  lattice  spacing  goes  to  zero  yields  the  following  analytical  approximation  of 
(22)-(24): 


xk 

Xk  /  /  /\  dx  , 

(25) 

yk 

Yk  =  j  /  i\dy' i  and 

(26) 

(27) 
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R-domain  r-domain 


Fig.  8  Example  of  the  spatial  transformation  of  an  operating  domain  (black)  and  a  sampling  tra¬ 
jectory  (blue)  generated  from  a  spanning-tree  coverage  algorithm  (red).  Reprinted  from  [54]  with 
permission  from  Elsevier. 


The  transformation  (25)-(27)  exists  and  is  invertible  if  ofa),  o(yk)>  and  t(^) 
are  positive,  bounded,  and  continuous  functions.  Either  the  numerical  or  analytical 
form  of  the  transformation  presented  above  can  be  used  for  sampling  a  nonstationary 
spatiotemporal  field,  but  the  numerical  transformation  has  the  advantage  that  it  is 
valid  for  any  function  of  the  decorrelation  scales  (as  long  as  they  are  continuous 
and  positive);  however,  it  has  (relatively)  stringent  conditions  for  invertibility  [54]. 
The  analytical  transformation  is  always  invertible  when  the  decorrelation  scales  are 
positive,  but  it  is  only  valid  when  the  spatial  scales  are  decoupled. 


4.2  Multi-vehicle  Coverage  Control 

This  section  describes  a  multi-vehicle  coverage  algorithm  that  invokes  the  coor¬ 
dinate  transformation  from  the  previous  section  to  generate  nonuniform  sampling 
trajectories  [54].  We  first  implement  a  multi- vehicle  controller  that  uniformly  cov¬ 
ers  the  ^-domain.  To  ensure  that  the  7?-domain  is  uniformly  covered  with  no  gaps 
or  overlaps  between  sensor  measurements,  we  employ  without  loss  of  generality  an 
existing  coverage  algorithm  called  Spanning  Tree  Coverage  (STC)  [21].  Intuitively, 
the  algorithm  works  like  a  lawn-mowing  algorithm  along  a  closed  path.  Path  closure 
is  crucial  due  to  the  temporal  nature  of  the  field;  locations  are  revisited  to  ensure  that 
the  temporal  variation  of  the  field  is  captured.  As  an  example,  consider  the  spatial 
scale  function  with  a  solitary  contraction  of  the  decorrelation  scales  and  the  sam¬ 
pling  boundary  in  Figure  7.  Figure  8  shows  the  spanning  tree  and  generated  path  in 
the  ^-domain  (left)  and  the  r-domain  (right). 
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Fig.  9  Mapping  error  of  a  3D  nonstationary  spatiotemporal  field  using  (a)  uniform  sampling  and 
(b)  a  nonstationary  sampling  algorithm.  Reprinted  from  [54]  with  permission  from  Elsevier. 


We  seek  to  steer  a  set  of  mobile  vehicles  to  motion  along  the  path  generated  by 
the  STC  algorithm.  The  vehicles  maintain  equal  spacing  from  one  another  along  the 
path.  To  illustrate  the  algorithm,  we  performed  a  simulation  of  twelve  vehicles  with 
random  initial  conditions.  The  trajectories  in  the  7?-domain  are  transformed  back 
to  the  r-domain  by  numerically  inverting  the  coordinate  transformation.  In  the  r- 
domain,  vehicles  cluster  their  measurements  in  the  area  of  the  dip  to  ensure  that  the 
full  domain  is  covered,  even  in  areas  where  the  spatial  and  temporal  decorrelation 
scales  are  small.  Figure  9  shows  the  resulting  mapping  error  in  comparison  to  a 
uniform  sampling  algorithm.  There  is  a  dip  in  the  decorrelation  scales  centered  at 
*0,1  =  yo,i  =0  km,  and  t  —  12  hours.  Five  vehicles  travel  along  a  closed  path  in  the 
r-domain.  The  blue  curve  on  the  t  =  0  plane  shows  the  spatial  path  that  the  vehicles 
track.  The  average  mapping  error  is  0.271  using  uniform  coverage  (Figure  7(a))  and 
0.178  when  using  the  nonstationary  sampling  algorithm  (Figure  7(b)).  Likewise,  the 
maximum  error  is  reduced  from  0.983  to  0.768. 


5  Summary 

This  report  describes  a  DDDAS  project  that  optimizes  vehicle  routing  using  two 
different  strategies.  The  first  strategy  is  based  on  maximizing  observability  of  an 
unknown  process,  e.g.,  a  flow  field.  This  strategy  is  applied  to  an  aircraft  refuel- 
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ing  scenario  in  which  the  trailing  aircraft  uses  distributed  pressure  sensors  to  center 
its  position  on  the  wake  of  the  lead  aircraft.  The  observability  of  the  leader  air¬ 
craft  is  assessed  using  measures  of  the  empirical  observability  gramian  and  a  re¬ 
cursive  Bayesian  filter  is  implemented  to  estimate  the  leader’s  wake  parameters.  A 
wavefront-propagation  optimal  control  algorithms  invokes  the  unobservability  in¬ 
dex  as  a  cost  metric  and  incorporate  estimates  of  the  leader’s  state  provided  by  a 
Bayesian  filter.  Numerical  simulations  of  formation  flight  and  aerial  refueling  ap¬ 
plications  illustrate  that  the  control  algorithm  successfully  steers  the  vehicle  to  a 
desired  relative  position  using  the  estimated  the  wake  parameters. 

The  second  strategy  for  DDDAS  vehicle  routing  is  based  on  a  statistical  repre¬ 
sentation  of  an  unknown  spatiotemporal  process.  This  strategy  is  applied  to  achieve 
non-uniform  coverage  of  a  spatiotemporal  environmental  field.  We  describe  a  multi¬ 
vehicle  control  algorithm  that  generate  nonuniform  coverage  of  a  nonstationary  spa¬ 
tiotemporal  field  by  using  a  coordinate  transformation  to  render  the  field  locally  sta¬ 
tionary  in  a  new  set  of  coordinates  (the  7?-domain).  Existing  spatial  coverage  and 
path  sampling  algorithms  produce  trajectories  that  uniformly  sample  the  /^-domain; 
in  the  original  coordinates  (the  r-domain),  these  trajectories  concentrate  measure¬ 
ments  near  regions  of  high  variability. 

The  DDDAS  concept  manifests  itself  in  both  applications  described  in  this  re¬ 
port  in  the  following  sense:  measurement  data  from  a  distributed  sensor  array  is 
assimilated  into  a  model  of  an  estimated  process,  which  is  in  turn  used  to  guide  the 
collection  of  subsequent  measurements.  The  report  describes  new  capabilities  with 
regards  to  observability-based  sensor  path-planning  and  adaptive  sampling  of  non¬ 
stationary  processes  for  which  non-uniform  coverage  is  required.  Ongoing  direc¬ 
tions  of  research  include  solving  the  co-array  problem:  i.e.,  simultaneous  mapping 
and  estimation  decorrelation  scales,  also  known  as  hyperparameters. 
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